/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined(_ParticlesOperators_)
#define _ParticlesOperators_
#include <src/miscellaneous/OS_basic.h>
#include <src/math/norms.h>
#include <src/math/Operators/weight_set_operator.h>

/*! \brief interaction for class ParticlesOperators
	This is the multidimensional case (dim >=2)
*/
template <class A, std::size_t Dim = A::Dim>
struct Particle_Interaction {
	typedef typename A::value_type value_type;
	typedef typename A::const_reference const_reference;
	typedef typename A::reference reference;
	//! fill an already initialized reference with \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. 
	void operator()(const_reference left, const_reference right, reference out) const
	{
		out = left;
		out -= right;
		out /= jmmath::norm22(out);
	};
	//! return \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. slow method, use 
	value_type operator ()(const_reference left, const_reference right) const
	{
		value_type a(left-right);
		return a/jmmath::norm22(a);
	};
};
/*! \brief interaction for one dimensional class ParticlesOperators
	This is the one dimensional case specialization of ParticlesOperators (Dim =1)
*/
template <class A>
struct Particle_Interaction<A,1> {
	typedef typename A::value_type value_type;
	typedef typename A::const_reference const_reference;
	typedef typename A::reference reference;
	//! fill an already initialized reference with \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. 
	void operator()(const_reference left, const_reference right, reference out) const
	{
		out = 1./(left-right);
	};
	//! return \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. slow method, use 
	value_type operator ()(const_reference left, const_reference right) const
	{
		return 1./(left-right);
	};
};

template <class T>
struct ParticlesOperators_interface 
{
	ParticlesOperators_interface() {};
	typedef typename T interfaced_type;
};

/*! \brief Base class for Particle operators

	Particle Operators are weighted_operator for which the interaction between two particles \f$ (S_i,S_j) \f$
	is \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. Thus the global behavior of the operator is
	\f[
		f( \{S_i \}_i = \{ \sum_{j \in \nu_i} w_{i,j}\frac{S_i-S_j}{|S_i-S_j|^2} \}_i
	\f].
*/
template <class A, class R >
class ParticlesOperators  : 
	public weighted_operator<A,R,Particle_Interaction<A> >
	, public virtual ParticlesOperators_interface<weighted_operator<A,R,Particle_Interaction<A> > >
{
public :
	typedef ParticlesOperators_interface												interface_type;
	typedef ParticlesOperators															type;
	typedef Particle_Operator															Family_Operator;
	typedef weighted_operator<A,R,Particle_Interaction<A> >								Base;
	typedef typename Base																Base_Cast;

	ParticlesOperators(){};
	ParticlesOperators(const weight_set_type & data) :  graph_base(data) {};
	ParticlesOperators(const Base &first) : Base(first)
	{};
	ParticlesOperators(const ParticlesOperators & from) :  Base(from) {
	};
	ParticlesOperators & operator = (const ParticlesOperators & from) {
		Base::operator =(from);
		return *this;
	};
	virtual ~ParticlesOperators(){};
	virtual smart_ptr_type Create_Instance() const {
		return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
	};
};

#endif